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c^ \ Abstract 

Partition of unity methods, such as the extended finite element method 
(XFEM) allow discontinuities to be simulated independently of the mesh [1]. 
^ , This eliminates the need for the mesh to be aligned with the discontinuity or 

CN ! cumbersome re-meshing, as the discontinuity evolves. However, to compute 

the stiffness matrix of the elements intersected by the discontinuity, a subdi- 
tJ- ■ vision of the elements into quadrature subcells aligned with the discontinuity 

t^ ■ is commonly adopted. In this paper, we use a simple integration technique, 

proposed for polygonal domains |2] to suppress the need for element sub- 
division. Numerical results presented for a few benchmark problems in the 
context of linear elastic fracture mechanics and a multi-material problem, 
show that the proposed method yields accurate results. Owing to its sim- 
plicity, the proposed integration technique can be easily integrated in any 
existing code. 
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1. INTRODUCTION 

The classical finite element method (FEM) is one of the clear choices to 
solve problems in engineering and science. But the classical FEM approaches 
fail or are computationally expensive for some classes of problems, viz., equa- 
tions with rough coefficients and discontinuities (arising e.g. in the reaction- 
diffusion equation, the advection-diffusion equation, crack growth problems, 
composites, materials with stiffeners etc.) and problems with highly oscilla- 
tory solutions viz., solution of Helmholtz's equation. In an effort to improve 
the FEM, Babuska et al., [3] showed that the choice of non-polynomial ansatz 
functions, when tailored to the problem formulation, lead to optimal conver- 
gence, whereas the classical FEM, relying on the approximation properties 
of polynomials, performs poorly. Also in the case of Helmholtz's equation, 
Melenk P] showed that plane waves displaying the same oscillatory behavior 
as the solution can serve as effective enrichment functions. This lead to the 
birth of the Partition of Unity Method (PUM). 

Belytschko's group in 1999 [H [S], exploited the idea of partition of unity 
enrichment of finite elements (Babuska et al., [3]) to solve linear elastic frac- 
ture mechanics problems with minimal remeshing. The resulting method, 
known as XFEM is classified as one of the partition of unity methods in 
which the main idea is to extend a classical approximate solution basis by 
a set of enrichment functions that carry information about the character of 
the solution (e.g., singularity, discontinuity, boundary layer). 

As it permits arbitrary functions to be incorporated in the FEM or the 
mesh-free approximation, partition of unity enrichment |6l [7] leads to greater 
fiexibility in modeling moving boundary problems, without changing the un- 
derlying mesh while the set of enrichment functions evolve (and/or their 
supports) with the interface geometry. 

In PU type methods, the enrichment is extrinsic and resolved through 
additional degrees of freedom. The enrichment can also be intrinsic, based 
on the recent work by Fries and Belytschko [8]. In this paper, we focus on 
the extrinsic partition of unity enrichment and in general, the field variables 
are approximated by[IlllllSl[iniEl[IIl[71[I2]: 

u (x) = 2_, ^/(x)q/ + enrichment functions (1) 



where Nj{x.) are standard finite element shape functions, q/ are nodal vari- 
ables associated with node I. 

XFEM, one of the aforementioned partition of unity methods, was suc- 
cessfully applied for crack propagation and other fields in computational 
physics [I31[T11|151II61[I71[IH1|191|201[2I1[22] and recently open source 
XFEM codes were released to help the development of the method [23] and 
numerical implementation and efficiency aspects were studied [21]. XFEM 
is quite a robust and popular method which is now used for industrial prob- 
lems [25] and under implementation by leading computational software com- 
panies. It is not the scope of this paper to review recent advances of partition 
of unity methods, and the interested readers are referred to the literature, 
for instance, Bordas and Legay [2S], Karihaloo et al., [27j and Belytschko et 
al, [28]. 

Although XFEM is robust and applied to a wide variety of moving bound- 
ary problems, the flexibility provided by this class of methods also leads to 
associated difficulties: 

• when the approximation is discontinuous or non-polynomial in an ele- 
ment, special care must be taken for numerical integration; 

• the low order of continuity of the solution leads to poor accuracy (esp. 
in 3D) of the derivatives close to regions of high gradient, such as crack 
fronts [29j which motivated recent work on adaptivity for GFEM [30| 
[31], meshfree methods [32l|33l[3l] and XFEM (Bordas and Duflot [35], 
Bordas et al. [36], Rodenas et al. [37]). 

An important flrst attempt to simplify numerical integration was by Ven- 
tura [38j, who focuses on the elimination of quadrature subcells commonly 
employed to integrate strongly or weakly discontinuous and non-polynomial 
functions present in the enriched FE approximation. His work is based on re- 
placing non-polynomial functions by 'equivalent' polynomials. The proposed 
method is exact for triangular and tetrahedral elements, but for quadrilateral 
elements, when the opposite sides are not parallel, additional approximation 
is introduced. 

Another method that alleviates this difficulty is strain smoothing [39|ll0]. 
The main idea is to combine the smoothed finite element method (SFEM) |4H 
W2\ with the XFEM. The SFEM relies on strain smoothing, which was pro- 
posed by Chen et al., [13] for meshless methods, where the strain is written as 
the divergence of a spatial average of the standard (compatible) strain field 



- i.e., symmetric gradient of the displacement field. In strain smoothing, the 
surface integration is transformed into an equivalent boundary integration by 
use of the Green-Ostrogradsky theorem. The Smoothed XFEM was intro- 
duced in [39], but much remains to be understood regarding the convergence, 
stability and accuracy of this method. 

In this paper, we propose to use the new numerical integration tech- 
nique proposed by the authors for arbitrary polygonal domains [21 SH SS] 
to compute the stiffness matrix. Each part of the elements that are cut or 
intersected by a discontinuity is conformally mapped onto a unit disk using 
Schwarz-Christoffel mapping. A midpoint quadrature is used to obtain the 
integration points as opposed to the regular GauB cubature rule. Thus, the 
proposed method which works only in 2D, eliminates the need to sub-divide 
the elements cut by discontinuities into quadrature subcells for the purpose of 
numerical integration of the stiffness matrix. The Schwarz-Christoffel map- 
ping has been applied to mesh free Galerkin method by Balachandran et 
al., [46j to obtain the weight function of the arbitrary shaped support do- 
main obtained from natural neighbor algorithm. 

The paper is organized as follows. In the next section, we briefly recall 
the basic equations of the XFEM. Section |3] will explain the new numerical 
integration scheme. The efficiency and convergence properties of the pro- 
posed method are illustrated in section H] with a few benchmark problems 
taken from linear elastic fracture mechanics and a multi-material problem, 
which is followed by some concluding remarks in the last section. 

2. BASICS OF PARTITION OF UNITY METHODS FOR DIS- 
CONTINUITIES AND SINGULARITIES 

In this section, we give a brief overview of partition of unity methods 
in FEM for problems with strong and weak discontinuities We focus on the 
extended FEM [H |2S], but the method is identical for alternatives such as 
theGFEM [301 [31]. 

With a regular finite element method, the mesh has to conform to the dis- 
continuities and a very fine mesh is required in regions with sharp gradients. 
When the discontinuity surface evolves, cumbersome remeshing is required. 
The XFEM alleviates these difficulties by allowing the discontinuities to be 
independent of the mesh. An XFEM model consists of a regular FE mesh, 
which is independent of the discontinuity geometry. Figure (JT]) illustrates a 
typical FE mesh with a crack. 







HH Ba enriched element, (tip element) 

^^ Reproducing elements 

H enriched element (split element) 

Partially enriched element (blending element) 

Standard element 

Figure 1: A typical FE mesh with a discontinuity. The squared nodes are enriched with 
the Heaviside function, and the circled nodes are enriched with near-tip asymptotic field 
obtained from Westergaard solution ^j 



The main idea is to extend the approximation basis by a set of enrichment 
functions, that are chosen based on the local behavior of the problem. For 
the case of linear elastic fracture mechanics, two sets of functions are used: 
a Heaviside jump function to capture the jump across the crack faces and 
asymptotic branch functions that span the two-dimensional asymptotic crack 
tip fields. The enriched approximation for fracture mechanics problems takes 
the form [H [201 HZ]: 



u'^(x)= J2 iV/(x)q/+5^iV^(x)iJ(x)aj+5^ iV^(x)5^i?,(x)b^, (2) 



a=l 



where aj and h^ are nodal degrees of freedom corresponding to the Heaviside 
function H and the near-tip functions, {-BQ}i<a<4, given by: 



{Ba}{r,6)i<a<4 = y/rl 
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(3) 

Nodes in set Jvf^ are such that their support is split by the crack and nodes in 
set 'N^ belong to the elements that contain a crack tip. These nodes are en- 
riched with the Heaviside and near-tip (branch functions) fields, respectively. 
In the discretization of Equation (I2]), the displacement field is global, but the 
supports of the enriching functions are local because they are multiplied by 
the nodal shape functions. 

This modification of the displacement approximation does not introduce 
a new form of the discretized finite element equilibrium equation, but leads 
to an enlarged problem to solve: 
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where the element stiffness matrix is given by: 
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where B^f^ is the standard strain-displacement matrix, Bg„^ and Bg„^ are 
the enriched parts of the strain-displacement matrix corresponding to the 
Heaviside and asymptotic functions, respectively. D is the material matrix. 
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Figure 2: Integration in an element with a straight and kinked discontinuity: standard 
decomposition of an element for integration of a discontinuous weak form for XFEM. 
Gaufi points are introduced within each triangle to ensure proper integration of the dis- 
continuous displacement field. 



The numerical integration of the stiffness matrix in elements intersected 
by a discontinuity, be it a material interface (weak discontinuity) or a crack 
(strong discontinuity) is not trivial. The standard GauB quadrature can- 
not be applied in elements enriched by discontinuous terms, because the 
Gaufi quadrature imphcity assumes a polynomial approximation. This is 



circumvented by partitioning the elements into sub cells aligned to the dis- 
continuity surface, in which the integrands are continuous and differentiable 
(see Figure ([2])). Although, the generation of quadrature subcells does not 
alter the approximation properties, it inherently introduces a 'mesh' require- 
ment. The steps involved in this approach are: (1) Split the element into 
subcells with the subcells aligned to the discontinuity surface, usually the 
subcells are triangular (see Figure ([2])) and (2) numerical integration is per- 
formed with the integration points from triangular quadrature. The subcells 
must be aligned to the crack or interface and this is costly and less accu- 
rate if the discontinuity is curved. Similar attempts were made to improve 
the integration of discontinuities in meshfree methods [IHl SH EUl ED [S2]- To 
alleviate this difficulty, we propose to use the new numerical integration tech- 
nique proposed by the authors [2] for polygonal domains to integrate over the 
elements intersected by the discontinuity. The next section will briefly review 
the Schwarz-Christoffel conformal mapping (SCCM) and then discuss how to 
numerically integrate over elements intersected by discontinuities using the 
SCCM. 

3. SCHWARZ-CHRISTOFFEL CONFORMAL MAPPING 

Conformal mapping is extremely important in complex analysis and finds 
its application in many areas of physics and engineering. A conformal trans- 
formation or biholomorphic map is a transformation that preserves local 
angles. In other words, if Fi and r2 are two curves that intersect at an angle 
9z in the z— plane at a point p, then the images / (Fi) and / (F2) intersect 
at an angle 9^, = 6^ at q = f (p). A Schwarz-Christoffel mapping is a trans- 
formation of the complex plane that maps the upper half-plane conformally 
to a polygon. 

Definition 1. A Schwarz-Christoffel map is a function f of the complex 
variable that conformally maps a canonical domain in the z— plane (a half- 
plane, unit disk, rectangle, infinite strip) to a 'closed' polygon in the w— 
plane. 

Consider a polygon in the complex plane. The Riemann mapping theorem 
implies that there exists a bijective holomorphic mapping / from the upper 
half plane {C G C : Im C > 0} to the interior of the polygon. The function / 
maps the real axis to the edges of the polygon. If the polygon has interior 
angles a, /3, 7, . . ., then this mapping is given by, 
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K 



dw 



(6) 



where i^ is a constant, and a < b < c < ... are the values, along the real 
axis of the ( plane, of points corresponding to the vertices of the polygon 
in the z plane. A transformation of this form is called a Schwarz-Christoffel 
mapping. 

Recently, the authors proposed a new numerical integration technique [21 
m US] using the Schwarz-Christoffel mapping and cubature rules on a disk 
to numerically integrate over arbitrary polygons that arise in polygonal finite 
element methods. Figure ([3]) shows the conformal mapping of an arbitrary 
polygon onto a unit disk on which a midpoint rule [53] or GauB Chebyshev 
rule [M] is used to obtain integration points. The distributions of integration 
points of the mid point quadrature and GauB- Chebyshev quadrature are 
illustrated in Figure (jl]). 



Conformal center 
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Figure 3: Mapping of the physical domain to the unit disk. This figure was produced with 
the MATLAB SC Toolbox M\ 




(a) Midpoint rule 




(b) Gaufi-Chebyshev rule 
Figure 4: Quadrature rules on a disk 
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In this paper, we propose to use conformal mapping in the context of 
the XFEM for 2D problems to numerically integrate over elements where the 
approximation or its derivatives is discontinuous. Each part of the element 
is conformally mapped onto a unit disk using the technique proposed in [2]. 
Figure ([5]) illustrates the above idea for split and tip elements. 



(a) 






SC mappin; 




Figure 5: Integration over an element with discontinuity (dotted line): (a) with kinked 
discontinuity, representing the split element and (b) strong discontinuity, representing the 
tip element. In both cases, the sub-polygon is mapped conformally onto the unit disk 
using Schwarz-ChristofFel conformal mapping. 

One of the cubature rules mentioned above is used to obtain integration 
points. For elements that are not intersected by a discontinuity surface, the 
standard isoparametric mapping is implemented with four integration points 
for each element. 
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4. NUMERICAL EXAMPLES 

In this section, we illustrate the effectiveness of the proposed method by 
solving a few benchmark problems taken from linear elastic fracture mechan- 
ics and a multi-material problem. We first consider an infinite plate under 
tension, then a plate with an edge crack, for which analytical solutions are 
available. Then, a multiple crack problem, followed by a multi-material prob- 
lem is considered and as a last example crack growth in a double cantilever 
beam is studied. In this study, unless otherwise mentioned, quadrilateral ele- 
mentqj are used. In case of the XFEM with a standard integration approach, 
13 integration points per subcell are used and a similar number of integration 
points are used for the SCCM, i.e., for a tip element with six subcells, 78 
(=13 x6) integration points are used for both methods. In this study, we 
use only topological enrichment, i.e., only the tip element is enriched by near 
tip functions [561 E3 l58] . 

4.I. Infinite plate under tension 

Consider an infinite plate containing a straight crack of length a and 
loaded by a remote uniform stress field a as shown in Figure ([6]). Along 
ABCD the closed form solution in terms of polar coordinates in a reference 
frame (r, 6) centered at the crack tip is 



^x{r,9) = —=cos-\l-sm-smYJ (7a) 

ay{r, ^) = — cos - ( 1 + sm - sm — 1 (7b) 

, ,. Kj . e 6 36 

Cxvir., o) = —= sm - cos - cos — (7c) 



The closed form near-tip displacement field is: 

u,{r,e) = A^U-^yfcos- 2 - 2z/ - cos^ - (8a) 

uy{r,9) = ^^§V^sin^ (2 - 2. - cos^ Q (8b) 



"Bilinear element, 4 noded quadrilateral element 
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Figure 6: Infinite cracked plate under remote tension: geometry and loads 



In the two previous expression Kj = a^/Tia denotes the stress intensity 
factor (SIF), v is Poisson's ratio and E is Young's modulus. All simulations 
are performed with a = 100mm and a = 10^ N/mm^ on a square mesh with 
sides of length 10mm. 

Before carrying out the mesh convergence study and other numerical 
studies, the influence of the number of integration points in the tip element 
on the numerical SIF is studied. A structured quadrilateral mesh (60 x 
60) is used for the study and the number of integration points are varied 
until the difference between two consecutive computations are less than a 
specified tolerance. During this study, the number of integration points for 
both methods are kept the same. The convergence of the SIF with the 
integration points is shown in Figure ([7]). It is seen that with the increase 
in the number of integration points, the SIF initially increases but reaches a 
constant value beyond 60 integration points. 

The convergence and rate of convergence in numerical stress intensity 
factor are shown in Figure ([H]). It is seen that for the same number of in- 
tegration points, the proposed method outperforms (although only slightly) 
the conventional XFEM. But with increase in mesh size, both techniques of 
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Figure 7: Griffith Problem: the convergence of the numerical stress intensity factor with 
number of integration points in the tip element. A structured quadrilateral mesh (60 x 60) 
is used. 



numerical integration approach the analytical solution. 

Examining Figure ([8]) shows that the convergence rates in the SIFs are 
suboptimal both for the XFEM with standard integration and the XFEM 
with the new integration technique proposed. This is for two reasons: (i) 
only the tip element is enriched (topological enrichment [2^ \5E[ \57\ , which 
asymptotically reduces the XFEM approximation space to the standard FEM 
approximation space, this limits the optimal convergence rate to 0.5 in the 
presence of a square root singularity; (ii)no blending correction [59] is per- 
formed, which leads to yet a smaller convergence rate of 0.4, which is consis- 
tent with the literature ^7\ ED] • 

To demonstrate the effectiveness of the proposed method in case of skewed 
elements, the domain is meshed with irregular elements. The coordinates of 
interior nodes are given by 



x' = X + (2rc — l)airAx 
y' = y+ (2rc - l)airAy 



(9a) 
(9b) 



where r^ is a random number between and 1.0, air is an irregularity factor 
controlling the shapes of the distorted elements and Ax, Ay are initial regular 
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Figure 8: Griffith Problem: the convergence of the numerical stress intensity factor to the 
analytical stress intensity factor and convergence rate. 
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element sizes in the x— and y— directions respectively. The discretization is 
shown in Figure (jH]). 







Figure 9: Polygonal mesh for infinite plate problem with a crack under uniform far field 
tension. The irregular mesh is generated with an irregularity factor, air = 0.4. 

The convergence of the numerical stress intensity factor for the irregular 
mesh is shown in Figure f lTOj) . However, the convergence for distorted ele- 
ments exhibits a non-uniform behavior. Results show that the results from 
both the numerical techniques are comparable. The advantage of the pro- 
posed method is that it eliminates the need to subdivide the split and the 
tip elements. 

4-2. Edge crack under tension 

A plate of dimension 1 x 2 is loaded by a tension a = 1 over the top edge. 
The displacement along the y-axis is fixed at the bottom right corner and the 
plate is clamped at the bottom left corner. The geometry, loading, boundary 
conditions and domain discretization are shown in Figure [TTl The reference 
mode I SIF is given by 



Ki 



avvra 



(10) 



where a is the crack length, H is the plate width and F{jj) is an empirical 
function given as (For {jj) < 0.6) 



-'i 



1.12-0.231 (-^) + 10.55 ('-^j -21.72 
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Figure 10: Griffith Problem: the convergence of the numerical stress intensity factor to 
the analytical stress intensity factor for distorted elements. 



The convergence of the mode I SIF with mesh size and the rate of con- 
vergence of the SIF for a plate with an edge crack is shown in Figure flT2|) . 
It is seen that with decrease in mesh density, both methods approach the 
analytical solution. Also, the proposed method performs slightly better than 
the standard XFEM. The rate of convergence for both methods is 0.4, very 
similar to the results available in the literature [U EOl [56] . The rate of con- 
vergence can be improved by using a 'geometrical' enrichment as suggested 
by Bechet et al., [56] or by modifying the enrichment functions such that 
they are zero in the standard elements and vary continuously in the blending 
elements as suggested by Fries (BT]. 

In all the above examples, the background FE mesh is made up of regular 
quadrilateral elements. The numerical integration of the stiffness matrix over 
regular quadrilateral elements is not a difficult task. The real challenge is 
when the background mesh is made up of polygons or when the crack faces are 
irregular, i.e., when the crack faces cut the elements in such a way that at least 
one of the subdomains, created by the intersection of the geometry with the 
mesh, is a polygonal with more than 3 edges. To demonstrate the usefulness 
of the proposed integration technique, we consider the problem of inclined 
crack in tension with bilinear quadrilateral element as background FE mesh. 
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Figure 11: Plate with edge crack under tension 
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Figure 12: Edge crack problem: the convergence of the numerical stress intensity factor 
to the analytical stress intensity factor and convergence rate. Both the methods show a 
rate of convergence of 0.4, very close to the ones obtained in [55] . 
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The study of problems involving polygonal FE meshes and arbitrary crack 
faces will be the topic of future papers. 

4-3. Inclined crack in tension 

Consider a plate with an angled crack subjected to a far field uniaxial 
stress field (see Figure (!T3|) ). In this example, Kj and Ku are obtained as 
a function of the crack angle /3. For the loads shown, the analytical stress 
intensity factors are given by [22] 

(12) 



i^7 = o"A/7racos/3cos/3, i^// = crv7racos/3sin/3. 
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Figure 13: Inclined crack in tension 

The influence of crack angle /3 on the SIFs is shown in Figure flT4|) . A 
structured mesh (100 x 100) is used for the study. For crack angles, < /9 < 
90, the crack face intersects the elements in such a way that one region of 
the split elements (either above or below the crack face) is a polygon. The 
polygonal subdomain is mapped onto a unit disk (see Figure ([5])) instead of 
subdividing it into triangles. It is seen from Figure ( !T^ that the numerical 
results are comparable with the analytical solution. 

4.4. Multiple cracks in tension 

In the next example, we consider a plate with two cracks. The geometry 
and boundary conditions of the problem are shown in Figure fITB]) . In this 
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Figure 14: Variation of stress intensity factors Ki and Ku with crack angle, /3. 

case, the problem is solved only by the new proposed method. The material 
properties are: Young's modulus £" = 3 x 10^ and Poisson's ratio, v = 0.3. A 
mesh size of 72 x 144 is used for the current study with crack size, 2ai = 0.2. 
The length of the other crack 2a2 is varied. 

Figure ( 1T6|) shows the variation of the mode I SIF for different H/L ratios 
and for different ratio of crack lengthcl with ^i = and ^2 = 0, where 9i 
and 02 are the angles subtended by the crack faces with the horizontal (see 
Figure (fT5|) ). The normalized mode I SIFu is plotted for point A in Figure 
( jT^ . This is done to non-dimensionalize the results. Also, the value of 
-^analytical IS taken as the value for a plate with a center crack given by 



Kr 



a\ / vrasec 



/ na 
\2^ 



(13) 



where a is the half crack length, w = ^ is the half width of the plate, and a 
is the far field tensile load apphed at the top of the plate. 

It can be seen that with increase in the H/L ratio, the interaction between 
the cracks decreases, which is as expected. It is also seen from Figure 
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Figure 15: Plate with multiple cracks: geometry and loads. The length of the cracks are 
2ai and 2a2. 9i and 62 are the angles subtended by the crack faces with the horizontal. 
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Figure 16: Influence of H/ L ratio and 2a2/2ai on computed stress intensity factor. 



( IT6|) that as the ratio of crack lengths increases, the normahzed mode I SIF 
also increases. And as the H/ L ratio increases, all the curves approach the 
analytical solution (-ft'anaiyticai)- The values will not be equal to the analytical 
solution, because the analytical solution is computed for the case with one 
center crack. 

Next, the influence of the relative angle between the cracks on the mode I 
and the mode II SIF for a crack length of 2ai = 2a2 = 0.2, with the distance 
between the cracks: H = 0.1 and L = 0.2 is studied. Figure ( 1T71) shows 
the variation of mode I and mode II SIF with the angle subtended by the 
cracks to the horizontal axis. It is seen that with the increase in the angle 
of the crack, {6i and 62), the mode I SIF decreases and approaches zero for 
61 = 62 = 90°. While the mode II SIF, initially increases with increase in the 
crack angle and reaches the maximum for the crack angle 61 = 62 = 45° and 
decreases with further increase in the angle. 



4-5. Bimaterial bar 

To illustrate the effectiveness of the proposed method to integrate weak 
discontinuities, we simulate a one- dimensional bimaterial bar discretized with 
2D quadrilateral elements. Consider a two-dimensional square domain Q = 
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(b) mode II 

Figure 17: Plate with two cracks: the variation of mode I SIF and mode II SIF with 
respect to angle between the cracks for crack lengths 2ai — 0.2 and 2a2 = 0.2. The 
distance between the cracks are: 7J = 0.1 and L = 0.2 
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r2iUf22 of length L = 2 with the material interface F located at 6 = L/2. The 
Young's modulus and Poisson's ratio in Qi = (— 1,&) x (—1, 1) are Ei = 1, 
z/ = 0, and that in Q.2 = {b, 1) x (—1, 1) are E2 = 10, z/ = 0. With no body 
forces, the exact displacement solution with Uy = at y = —1 and My = 1 at 
y = 1 is given by: 

, (y + l)a, -l<y<b, 

u{y)={ ^^ / ' -'- ' (14) 



1 + f (y - 1) a, b<y<l 



where, 



« = z..., .^ z... .^ (15) 



E2{b+l)-Ei{b 
The relative error in displacement norm used to measure the accuracy of 



the results: 
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Figure (ITSl) shows the relative error and the rate of convergence in the 
displacement norm. It is seen that with decrease in mesh density, the relative 
error in displacement norm also decreases. The rate of convergence is also 
shown in the Figure (1181) . The proposed method obtains convergence rate of 
1.34 in the L2— norm, very close to the value reported in the literature [6T] . 
The convergence rate is suboptimal due to the absence of treatment of the 
blending elements [591I291I561I57] • 

Next, the influence of the material interface is studied. Three different 
configurations of the material interface are considered for the study (see 
Figure (^Uij). 

Figure (120]) shows the strain energy convergence with mesh refinement 
for three different configurations. It is seen that the strain energy converges 
as the mesh is refined. 

4-6. Double cantilever beam 

The dimensions of the double cantilever beam (see Figure fl?Tl) ) are 6x2 
and the initial pre-crack with length of a = 2.05 is considered. The material 
properties are taken to be Young's modulus, E = 100 and Poisson's ratio. 
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Figure 18: Bimaterial problem: convergence in the displacement (L2) norm: (a) the 
relative error and (b) the convergence rate. The method obtains convergence rate of 1.34 
in the L2— norm 
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Figure 19: Bimaterial probleni: (a)Straight Interface; (b) Slanted Interface (Positive slope) 
and (c) Slanted Interface (Negative slope) 
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Figure 20: Bimaterail problem: convergence of the strain energy with mesh refinement 

V = 0.3. And the load P is taken to be unity. By symmetry, a crack 
on the mid-plane of the beam is under pure mode I and the crack would 
propagate in a straight line, however, due to small perturbations in the crack 
geometry, the crack takes a curvilinear path p. A quasi-static crack growth 
is considered in this study and the growth is governed by the maximum hoop 
stress criterion [63], which states that the crack will propagate from its tip 
in the direction 6c where the circumferential (hoop) stress aeg is maximum. 
The critical angle is computed by solving the following equation: 



Kj smiOc) + Kjj{3cos{ec) - 1) = 
Solving Equation (TT7|) gives the crack propagation angle 



(17) 



6r = 2 arctan 



-2(ff 



1 






(18) 



The crack growth increment, Aa is taken to be 0.15 for this study and the 
crack growth is simulated for 8 steps. The domain is discretized with a 
structured mesh consisting of 1200 elements. The crack path is simulated 
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t p = l 



a = 2.05 



D = 2 









L = 6 

I ^ 

P = l 

Figure 21: Geometry and loads of a double cantilever beam 



The crack path quahtatively 



using both methods and is shown in Figure (122 
agrees with the pubhshed results [1]. 

5. CONCLUSION 



In this paper, we used the new numerical integration proposed for arbi- 
trary polygons in [2] to integrate the discontinuous and singular integrands 
appearing in the XFEM stiffness matrix. The proposed method eliminates 
the need to sub-divide elements cut by strong or weak discontinuities or con- 
taining the crack tip. With a few examples from linear elastic fracture me- 
chanics and a bimaterial problem, the effectiveness of the proposed method is 
illustrated. It is seen that for similar number of integration points, the pro- 
posed technique slightly outperforms the conventional integration method 
based on sub-division. With mesh refinement, both integration techniques 
provide convergence of the SIFs to the analytical SIFs. It seems possible 
that the proposed technique could serve as a way to integrate discontinuous 
approximations in the context of 3D problems as well, which will be the topic 
of forthcoming communications. 
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Figure 22: Double cantilever beam: comparison of crack path between the two numerical 
integration methods. 
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